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ABSTRACT 


This  report  discusses  the  extension  of  an  infiltration  predicting  technique  to 
the  prediction  of  inter-room  air  movements.  The  air  flow  through  openings  is 
computed  from  the  ASHRAE  crack  method  together  with  a mass  balance  in  each 
room.  Simultaneous  solution  of  the  mass  balances  in  all  rooms  having  both 
large  and  small  openings  is  accomplished  by  a slightly  modified  Newton's  method. 

A simple  theory  for  two-way  flow  through  large  openings  is  developed  from  con- 
sideration of  density  differences  caused  by  different  temperatures  in  adjoining 
rooms.  The  technique  is  verified  by  comparison  to  published  experimental 
results.  The  results  indicate  that  the  simple  model  provides  reasonable  results 
for  complex  two  way  flows  through  openings.  The  model  is  as  accurate  as  the 
available  data,  that  is,  about  +20%.  The  air  flow  algorithm  allows  infiltration 
and  forced  air  flows  to  interact  with  the  doorway  flows  to  provide  a more  general 
simulation  capability. 
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building  energy  analysis;  building  heat  transfer;  computer  modeling; 
convection;  infiltration;  ventilation. 


PREFACE 


This  report  is  one  of  a series  documenting  MBS  research  and  analysis  efforts  in 
developing  energy  and  cost  data  to  support  the  Department  of  Energy/National 
Bureau  of  Standards  Measurements  Program.  It  was  prepared  by  the  Thermal  Anal- 
ysis Group,  Building  Physics  Division,  Center  for  Building  Technology,  National 
Engineering  Laboratory,  National  Bureau  of  Standards  (NBS).  This  work  was 
jointly  sponsored  by  (NBS).  The  development  of  multi-room  air  flow  modeling 
was  supported  by  DoE/NBS  Task  Order  A008  under  Interagency  Agreement  No.  E-77- 
A-01-6010.  This  report  describes  a computer  program  which  was  written  as  part 
of  an  effort  to  develop  a comprehensive  modeling  technique  for  predicting  the 
simultaneous  transfer  of  air,  moisture,  and  heat  in  and  through  multi-room 
buildings . 

The  author  wishes  to  acknowledge  the  support  and  direction  given  by 
Dr.  T.  Kusuda,  Thermal  Analysis  Group  Leader. 
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1 .  INTRODUCTION 


1 . 1 BACKGROUND 


Although  numerous  building  thermal  modeling  techniques  and  computer  programs , 
for  example  NBSLD  [1],  BLAST  [2],  and  DOE-2  [3],  exist  throughout  the  United 
States,  none  of  the  existing  techniques  /'programs  handle  the  following  processes 
simultaneously: 

- envelope  heat  transfer 

- envelope  air  leakage 

- envelope  solar  heat  gain 

- room-to-room  heat  transfer 

- room-to-room  air  and  moisture  transfer 

- intra-room  air  movement 

- energy  consumption  by  the  heating/cooling  equipment 

- indoor  comfort 

- water  vapor  condensation  and  contaminant  migration. 

Existing  models  are  virtually  all  single-room  models  where  dynamic  coupling 
between  the  heated  and  non-heated  spaces  and/or  the  cooled  and  non-cooled 
spaces  are  ignored. 

Comprehensive  multi-room  building  simulation  capabilities  will  be  needed  in 
the  coming  years  for  the  following  reasons: 

1.  Intra-room  convection  plays  a significant  role  not  only  for  the  transfer 

of  heat  from  the  interior  surfaces  to  the  room  air  but  also  for  the  thermal 
comfort  of  the  occupants.  Yet  existing  computational  technology  for  pre- 
dicting the  room  temperature  stratification  and  room  air  motion  is  very 
inadequate . 

2.  Passive  solar  design  techniques  are  expected  to  be  used  to  a large  extent 
in  new  building  design.  Proper  design  is  going  to  depend  on  understanding 
and  being  able  to  predict  natural  energy  flows  within  the  buildings. 

3.  Proper  control  of  the  air  flow  in  the  central  air  system  requires  accurate 
knowledge  of  room-by-room  energy  demand. 

4.  Moisture  and  contaminant  distribution  throughout  a building  must  be  capable 
of  being  predicted  in  order  to  insure  adequate  designs  as  the  emphasis  on 
"tighter"  buildings  increases. 

5.  Until  such  capability  exists,  effective  design  of  internal  partitioning 
with  respect  to  natural  ventilation  can  only  be  done  by  empirical  techniques. 

In  response  to  these  needs,  a research-oriented  computer  program  has  been 
developed  to  allow  the  detailed  study  of  simultaneous  air,  heat,  and  moisture 
transfer  in  and  through  a building  with  complex  internal  architecture.  This 
program  is  called  the  Thermal  Analysis  Research  Program  (TARP).  Documentation 
[4]  for  the  program  will  soon  be  published.  Primary  emphasis  in  the  develooment 
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of  TARP  has  been  on  air  transfer  because  this  is  a major  shortcoming  of  present 
techniques  and  because  it  is  basic  to  further  developments  in  moisture  and 
contaminant  analysis.  A previous  report  [5]  describes  initial  results  in  the 
development  of  TARP  and  listed  five  significant  areas  for  further  research: 

1.  Calculation  of  air  flows  through  large  openings  in  reasonable 
computation  time, 

2.  Calculation  of  two-way  flows  between  rooms, 

3.  Accurate  evaluation  of  the  wind  induced  pressure  distribution  around 
the  envelope  of  the  building, 

4.  Calculation  of  the  effects  of  room  air  stratification,  and 

5.  Availability  of  data  for  estimating  the  opening  areas  in  the  envelopes 
of  buildings. 


This  report  will  address  the  first  two  areas. 
1.2  OBJECTIVE 


The  overall  objective  of  this  research  effort  is  to  develop  a comprehensive 
modeling  technique  for  predicting  simultaneous  transfer  of  air,  moisture,  and 
heat  in  and  through  multi-room  buildings.  The  objective  of  this  particular 
study  is  the  development  of  a method  for  predicting  air  flows  through  large 
openings  between  rooms. 

1.3  SCOPE 


Section  2 of  this  report  discusses  the  extension  of  an  infiltration  predicting 
technique  to  the  prediction  of  inter-room  air  movements.  The  technique  is  veri- 
fied by  comparison  to  published  experimental  results  in  section  3.  Conclusions 
and  recommendations  are  given  in  section  4. 
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2.  CALCULATION  METHODS 


2.1  THE  AIR  FLOW  EQUATION 


The  TARP  air  flow  algorithm  is  described  in  detail  in  reference  [5]  where  the 
program  was  called  the  Multi-Room  Loads  Program  (MRLP).  The  air  flow  algorithm 
is  based  on  the  equation  [6] 

F = K * (AP)X  (1) 


where 


F = flow  rate  (kg/s) 

K = a constant 

AP  = pressure  difference  across  an  opening  (Pa) 

X = flow  exponent 

Pressure  differences  arise  from  wind,  air  density  differences,  and  system 
induced  flows.  By  applying  equation  (1)  to  all  openings  in  the  building 
envelope  and  all  openings  between  rooms  and  requiring  a mass  balance  in  each 
room,  a set  of  simultaneous  non-linear  algebraic  equations  is  created  which 
can  be  solved  for  the  zone  pressures  and  the  air  flow  through  each  opening. 

An  estimate  for  the  value  of  K was  made  by  referring  to  the  orifice  equation: 

F=C*A*p*  / 2 * AP  / p (2) 


where 


C = flow  coefficient 
A = observed  opening  area  (nP) 
p = air  density  (kg/m^) 

When  the  opening  is  small,  C equals  0.6  for  a wide  range  of  Reynolds  numbers 
as  shown  in  figure  1 [7],  TARP  assumes  this  value  of  C as  a default.  TARP 
allows  a variable  flow  exponent  (instead  of  the  0.5  of  the  orifice  equation) 
because  building  presurization  measurements  typically  correlate  to  equation 
(1)  when  X equals  about  0.65  [5]. 

2.2  SOLUTION  OF  THE  FLOW  EQUATIONS 


The  development  of  a solution  technique  for  the  air  flow  equations  has  been 
particularly  troublesome.  Efforts  have  focused  on  two  techniques  which  are 
described  by  Conte  and  DeBoor  [8].  The  first  technique  is  the  classic  Newton's 
method.  The  mass  balance  requirement  may  be  expressed  as 

l V±  = 0 (3) 

for  every  room.  The  flows,  Fj_ , are  summed  over  all  openings,  i,  in  the  enclosing 
surfaces  of  the  room.  In  Newton's  method,  successive  values  of  room  pressure, 

Pn,  for  each  room,  n,  are  calculated  by 
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(4) 


(k+1)  (k)  _ 

where  k is  the  iteration  number  and  C is  computed  from  the  matrix  relationship 

I J ] * [ C ] - [ B ] (5) 

where  [ B ] is  a column  matrix  each  element  given  by 

Bn  = Z F±  (6) 

and  J is  the  square  (i.e.,  N by  N for  a building  of  N rooms)  Jacobian  matrix. 
The  values  of  the  diagonal  elements  of  the  Jacobian  matrix  are  given  by 


J 


n 


n 


Z 

i 


3Fi 

3Pn 


(7) 


for  all  openings,  i,  into  room  n.  The  values  of  the  other  elements  are  given 
by 


J 


n , m 


Z 

i 


(8) 


for  all  openings,  i,  between  room  m and  room  n.  Iteration  proceeds  until  the 
net  air  flow  into  every  room,  Bn,  is  sufficiently  close  to  zero.  These  partial 
derivatives  are  very  easy  to  compute : 


3F  • 

— - = X,  * F • / AP  (9) 

8Pn  1 1 


d¥± 

— - = X,  * F,/AP 

3Pm 

m 


or 


3F  • 


3P 


1 = 0 


m 


(inter-room  surface)  (10) 


(envelope  surface)  (11) 


Note  that  the  term  (AP)‘^  ^ has  been  eliminated  from  these  expressions,  thus 
allowing  evaluation  of  the  derivatives  by  a simple  division  rather  than  a time 
consuming  exponential.  Also  note  that  as  AP  approaches  zero,  the  derivative 
is  undefined. 


The  second  method  calls  for  successively  approximating  each  room  pressure 
according  to 
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(12) 


p (k+l)  = p (k)  _ zy  /T.  — - 

n n 1 3Pn 

where  k is  again  the  iteration  number.  This  method  is  quite  simple  and  requires 
less  storage  space  than  Newton’s  method  because  it  does  not  use  the  Jacobian 
matrix. 

Initial  tests  of  these  two  methods  showed  that,  although  Newton’s  method  was 
fastest  for  most  test  cases,  it  occasionally  would  not  converge.  The  second 
method  converged  for  the  original  test  cases  and  was  chosen  for  the  initial 
version  of  TARP  [5],  It  was  then  found  that  this  method  converged  very  slowly 
when  the  openings  between  rooms  were  much  larger  than  the  openings  in  the 
building  envelope  which  is  a very  common  condition.  This  problem  led  to  a re- 
examination of  the  Newton's  method.  A simple  four  room  test  case  was  found  to 
usually  be  quadratically  convergent  (about  4 iterations)  except  for  a few  cases 
where  it  converged  verly  slowly  (about  30  iterations).  In  those  cases  it  was 
found  that  successive  iterations  were  over-correcting.  That  is,  they  were  suc- 
cessively far  above  and  then  far  below  the  correct  solution  because  successive 
corrections  were  of  about  the  same  magnitude  but  opposite  sign.  Convergence 
could  be  achieved  by  reducing  the  size  of  the  pressure  correction  by  about  one 
half.  Since  the  Newton's  method  is  normally  rapidly  convergent,  it  is  also 
necessary  to  reduce  the  size  of  the  correction  only  when  over-correcting  occurs. 
The  reason  for  occasional  slow  convergence  of  the  Newton’s  method  has  not  been 
found.  It  can  occur  with  nothing  more  than  a change  in  wind  direction  from 
what  was  a quardratically  convergent  case.  It  often  occurs  when  the  wind  and 
stack  pressures  are  about  equal.  Convergence  is  always  fastest  when  the  flow 
exponent  for  inter-room  openings  is  near  one.  The  experimental  studies  below 
indicate,  however,  that  the  exponent  should  be  one  half. 

The  air  flow  algorithm  was  incorporated  into  a test  program  which  allowed 
various  solution  techniques  to  be  studied  without  revising  TARP,  which  is  very 
large.  A listing  of  the  test  program  is  attached  as  appendix  A.  Details  of 
the  solution  algorithm  are  discussed  in  that  listing. 

The  Newton's  method  requires  the  simultaneous  solution  of  matrix  equation  (5) 
at  each  iteration.  Several  techniques  were  considered  for  that  solution.  The 
first  choice,  and  the  one  ultimately  chosen,  was  Gauss  elimination.  It  has  the 
disadvantage  that  solution  time  is  proportional  to  the  cube  of  the  number  of 
equations,  N,  when  N is  large.  The  number  of  equations  is  equal  to  the  number 
of  rooms.  The  Cholesky  method  [9]  is  somewhat  faster  at  large  N,  but  is  was 
found  to  be  sensitive  to  computer  truncation  errors.  A Gauss-Seidel  iteration 
was  also  tried  since  it  has  solution  time  proportional  to  the  square  of  N. 
However,  this  iteration  also  failed  for  large  openings  in  much  the  same  way 
that  equation  (12)  did.  The  Newton's  method  was  tested  for  larger  numbers  of 
rooms.  It  was  found  that  the  number  of  iterations  increased  with  the  number  of 
rooms.  Thus,  many  room  require  both  longer  iterations  and  more  iterations. 

The  number  of  iterations  was  also  found  to  increase  with  the  size  of  the  inter- 
room openings,  but  it  did  not  increase  as  dramatically  as  it  had  done  with 
equation  (12).  Therefore,  the  Newton's  method  is  most  appropriate  for  a small 
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number  of  rooms  with  large  execution  time  penalties  paid  for  simultaneously 
solving  many  rooms. 

TARP  uses  an  hourly  heat  balance  in  estimating  dynamic  room  energy  requirements. 
This  heat  balance  is  solved  iteratively,  and  at  each  iteration  a quasi-steady 
solution  of  the  air  flow  equations  is  obtained  by  the  Newton's  mehtod  described 
above.  Techniques  used  to  reduce  the  number  of  heat  balance  iterations 
contribute  to  the  overall  efficiency  of  TARP  and  are  described  in  [5]. 

2.3  A THEORY  OF  FLOW  THROUGH  LARGE  OPENINGS 

Equation  (1)  permits  air  to  flow  only  one  direction  through  an  opening.  Large 
openings,  such  as  doorways,  may  have  two  way  flow  as  the  stack  effect  between 
two  rooms  may  cause  a positive  AP  at  the  bottom  of  the  doorway  and  a negative 
AP  at  the  top  (or  vice  versa).  A theoretical  estimate  of  the  stack  induced  air 
flow  through  a large  opening  in  a vertical  partition  (a  doorway)  is  given  by 
Brown  and  Solvason  [9].  The  following  discussion  shows  that  the  TARP  method  is 
equivalent.  Figure  2 shows  a cross  section  of  a rectangular  opening  of  height 
H and  width  W in  a vertical  partition  which  separates  two  rooms  at  temperatures 
Tt  and  T^.  The  absolute  pressure  at  the  centerline  (z  = 0)  is  everywhere  equal. 
The  pressure  difference  caused  by  stack  effect  at  height  z is 

pl  " p2  = (Pi  - P2)  * g * z (13) 

The  volumetric  flow  through  an  infinitesimal  area  is  given  by 

dQ  = C * W * dz  * / 2 * Ap/p  (14) 

which  can  be  integrated  to  give  the  flow  through  the  top  half  of  the  opening 

Z=H/2  

Q = / dQ  = C/3  * W * / g * Ap/p  * H3  (15) 

Z=0 


The  coefficient  of  thermal  expansion  is  p = -Ap/(p  * AT)  where  p is  the  average 
density.  For  computational  simplicity  TARP  uses  the  density  of  the  incoming 
fluid  instead,  but  this  cannot  cause  a significant  error  at  normal  temperatures. 


Other  definitions  are: 

heat  transfer  rate:  q = Q*p*c*  (Tj_  - T^)  (16a) 
heat  transfer  coefficient:  h = q / [W  * H * (T^  - T2)]  (16b) 
Nusselt  number:  Nu  = h * H / k (I6c) 
Prandtl  number:  Pr  = c * u /k  (16d) 
Grashof  number:  Gr  = p“*g*|9*(T^-T2)*H^/u  (16e) 


The  expressions  can  be  substituted  into  equation  (15)  to  give 
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Nu  = C/3  * Pr  * / Gr 


(17) 


This  simplified  analysis  has  neglected  viscous  effects  and  effects  of  an  air 
velocity  parallel  to  the  partition.  According  to  Brown  and  Solvason  the  viscous 
effects  reduce  the  air  flow  through  thick  partitions  and  an  air  velocity  may 
either  increase  or  reduce  the  air  flow. 

TARP  can  handle  the  two-way  air  flow  through  a doorway  by  dividing  the  door 
into  an  upper  and  a lower  opening.  It  is  easy  to  determine  the  heights  of  the 
two  openings  which  cause  a stack  effect  giving  the  same  volumetric  flow  as  the 
Brown  and  Solvason  model.  These  turn  out  to  be  13/18  * H for  the  upper  opening 
and  5/18  * H for  the  lower  one. 

Brown  also  studied  openings  in  horizontal  partitions  [10].  He  developed  a 
theoretical  expression  for  convection  through  such  openings: 

Nu  = (0.29  to  0.35)  * Pr  * f~Gr  (18) 

In  this  case  the  Nusselt  and  Grashof  numbers  are  based  on  the  thickness  of  the 
partition.  This  thickness  is  the  vertical  space  available  for  the  development 
of  a stack  effect. 
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3.  VALIDATION  OF  THE  LARGE  OPENING  MODEL 


3.1  FLOW  THROUGH  OPENINGS  IN  VERTICAL  PARTITIONS 


Weber  and  Kearney  [11]  give  a correlation  for  the  two  way  flow  through  a doorway 
based  on  temperature  measurements  in  the  doorway 


Nu  = 0.26  * Pr  * / Gr 


(19) 


and  another  based  on  average  room  temperatures 


Nu  = 0.3  * Pr  * / Gr 


(20) 


These  correlations  are  based  on  similitude  experimental  studies  of  the  room 
shown  in  figure  3 and  tests  on  a full  scale  structure  with  a similar  configura- 
tion. Weber  and  Kearney  estimated  that  they  should  be  dependable  to  within  20 
percent.  They  also  compare  their  results  to  two  other  studies  (figure  4) 
including  Brown  and  Solvason's.  Comparing  equations  (19)  and  (17)  gives  a 
value  of  0.78  for  the  flow  coefficient,  C.  Equation  (20),  which  uses  average 
room  temperatures,  gives  a value  of  0.90  for  the  flow  coefficient.  This  seems 
unreasonably  large  and  is  apparently  due  to  the  uneven  temperature  distribu- 
tions which  occur  in  real  rooms,  especially  above  the  door  openings.  A study 
of  a dorrway  between  a small  room  and  a much  larger,  high  heat  loss,  two-story 
room  did  not  agree  well  with  the  correlation  (20). 

The  TARP  model  dividing  the  doorway  into  halves  was  used  for  a wide  range  of 
parameters  and  the  results  compared  to  equation  (19).  The  inter-room  mass  and 
energy  flows  were  computed  for  five  values  of  AT  (2,  4,  6,  8,  and  10  C)  each 
at  five  doorway  heights  (1.0,  1.5,  2.0,  2.5,  and  3.0  meters).  A flow  coeffi- 
cient of  .78,  a flow  exponent  of  0.5,  and  stack  heights  at  5/18ths  and  13/l8ths 
of  the  doorway  height  give  TARP  computed  mass  and  energy  flows  which  agree  with 
equation  (19)  to  within  0.1  percent  for  all  cases. 

3.2  FLOW  THROUGH  OPENINGS  IN  HORIZONTAL  PARTITIONS 

Figure  5 shows  the  experimental  results  (H  = thickness;  L = width  of  square 
opening)  of  Brown's  study  of  openings  in  horizontal  openings  [10].  In  this 
configuration  there  is  a significant  frictional  effect  in  the  thicker 
partitions.  Brown  incorporated  this  into  a single  equation: 


For  a TARP  evaluation  it  is  better  to  to  consider  the  thickness  effect  as  a 
modifier  to  the  stack  height.  A TARP  model  of  an  opening  in  a horizontal 
partition  would  again  divide  it  into  two  openings,  one  a distance  Z above  the 
center  of  the  partition  and  one  and  equal  distance  below  it.  Then,  for  C = 
0.78,  the  values  of  Z for  different  H/L  would  be: 


Nu  = 0.0546  * pr  * Gr-55  * (L/H)*33 


(21) 


8 


H/L 

Z 

0833 

. 168*H 

167 

. 120*H 

333 

. 093*H 

667 

.074*H 

These  almost  exactly  duplicate  the  Brown  results. 

There  are  still  several  questions  about  this  model.  The  effects  beyond  the 
studied  ranges  of  Gr  and  H/L  are  not  known.  The  model  predicts  that  the  flow 
should  go  to  zero  as  H approaches  zero!  There  should  be  some  flow.  It  is 
possible  that  two  separate  openings  would  behave  differently  than  a single 
opening  of  equivalent  area  because  one-way  flow  could  develop  in  each.  These 
questions  indicate  room  for  further  experimental  work.  In  addition,  a TARP 
model  would  have  to  allow  for  no  flow  between  rooms  when  the  upper  is  warmer 
than  the  lower. 
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4.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  results  indicate  that  the  simple  TARP  model  provides  reasonable  results  for 
complex  two  way  flows  through  openings.  The  model  is  as  accurate  as  the  avail- 
able data,  that  is,  about  +20  percent.  The  TARP  air  flow  algorithm  allows 
infiltration  and  forced  air  flows  to  interact  with  the  doorway  flows  to  provide 
a more  general  simulation  capability.  Although  more  work  must  be  done  for  other 
building  configurations,  the  close  match  between  the  TARP  algorithm  and  the 
experimental  correlations  is  encouraging. 

Because  of  the  rapid  increase  of  calculation  time  with  the  number  of  rooms 
simulated,  it  is  recommended  that  large  buildings  with  many  rooms  be  simulated 
by  dividing  the  building  into  groups  of  closely  coupled  rooms.  First,  treat 
each  group  of  rooms  as  a single  room  to  solve  for  the  infiltration  through  the 
building  envelope  and  air  flows  between  the  groups  of  rooms.  Then,  while  treat- 
ing those  air  flows  as  constant  gains  or  losses  to  the  appropriate  rooms,  solve 
for  the  air  flows  between  the  individual  rooms  in  each  group. 

Further  study  is  needed  on  the  simulation  of  room  air  stratification  both  for 
its  direct  effects  on  comfort  and  energy  requirements  as  well  as  for  its  effect 
on  inter-room  air  movement.  Study  is  also  needed  on  the  calculation  of  the 
wind  pressure  on  the  building  envelope.  Then  detailed  validation  should  be  per- 
formed with  carefully  selected  full  scale  tests.  After  successful  validation, 
it  will  still  be  necessary  to  develop  air  openings  data  for  typical  construction 
components  and  building  techniques  for  the  analysis  and  design  of  buildings. 
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Figure  1.  Orifice  flow  ueter  coefficients  [6] 
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Figure  .6.  Experimental  results  for  air  flow  through  a square 
opening  in  a horizontal  partition  [10] 


APPENDIX  A:  LISTING  OF  AIR  FLOW  TEST  PROGRAM 


The  subroutines  and  variables  used  in  this  program  were  prepared  for  easy 
inclusion  in  TARP . Several  features  could  have  been  simplified  or  improved  if 
the  only  goal  had  been  to  develop  a stand-alone  air  flow  analysis  computer 
program. 

Descriptions  of  variables  in  common: 


ACNVG1 

ACNVG2 

ACNVG3 

AFSPTR 

AIRDEN 

AMAXIT 

CPAIR 

DIR 

DTR 

FAHS 

FAREA 

FAZM 

FEXP 

FRACT 

HCOUNT 

LIST 

MAXAFS 

MAXZON 

MCPM 

MCPTM 

NAFS 

NZON 

OBP 

ODB 

PS 

PW 

PZ 

SPD 

SORTDZ 

STDTIM 

SUMAF 

TZ 

zs 

ZT 

ZZ 


relative  air  flow  convergence, 
minimum  air  flow  (kg/s), 
minimum  pressure  difference  (Pa). 

air  flow  opening  pointers  (near  side  and  far  side  room  numbers). 

ambient  air  density  (kg/m^). 

maximum  air  flow  mass  balance  iterations. 

specific  heat  of  air  (kWh/kg  C). 

wind  direction  (degrees). 

conversion  factor:  degrees  * DTR  = radians, 

air  flow  from  the  air  handling  system  (m^/s). 
effective  flow  area  (ra-). 
azimuth  angle  of  surface  (degrees), 
air  flow  exponent. 

fraction  of  computed  correction  (to  speed  convergence), 
count  of  heat  balance  iterations. 

printed  output  control  flag  (0  = least  detail,  2 = most). 

parameter:  maximum  number  of  air  flow  openings. 

parameter:  maximum  number  of  zones. 

inward  mass  flow  rate  time  specific  heat  (W/C) 

inward  energy  flow  of  moving  air  (W) . 

number  of  air  flow  openings. 

number  of  zones  (rooms). 

ambient  barometric  pressure  (Pa). 

ambient  dry  bulb  temperature  (C). 

stack  pressure  (Pa) . 

wind  pressure  (Pa). 

zone  (room)  air  pressures  (Pa).  [zone  0 = ambient] 
wind  speed  (m/s)  

square  root  of  the  zone  air  density  (/kg/m^)  . [zone  0 = ambient] 

hour  counter  (standard  time). 

sum  of  (unsigned)  zone  air  flows  (kg/s). 

zone  air  temperature  (C).  [zone  0 = ambient] 

height  of  opening  (m) . 

maximum  height  of  surface(m). 

height  of  zone  (m) . [zone  0 = ambient] 


Listing  of  subroutines: 


The  driver  program,  MAIN,  is  structured  to  test  convergence  for  various 
convergence  limits  and  ambient  conditions  for  a single  building  configuration. 
Other  driver  programs  were  written  to  test  air  flows  through  openings  for 
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comparison  to  the  experimental  results.  This  routine  uses  NAMELIST  input  for 
the  control  variables.  A sample  input  deck  is  shown  in  Appendix  B. 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 


C 


c 


c 


c 


PROGRAM  MAIN 

COMMON  /ENVT/  ODB( 1 ) ,0BP( 1 ) ,SPD( 1 ) ,DIR( 1 ) , AIRDEN, CPAIR, DTR, STDTIM, 
LIST, HCOUNT, ACNVG1 , ACNVG2 , ACNVG3 , AMAXIT , FRACT , NAFS , NZON 

INTEGER  STDTIM,  HCOUNT,  AMAXIT,  NAFS,  NZON 
REAL  ODB , OBP , SPD , DIR,  AIRDEN,  CPAIR,  DTR, 

ACNVG1,  ACNVG2 , ACNVG3 , FRACT 

NAMELIST  /CONTROL/  AMAXIT, ACNVG1 ,ACNVG2 ,ACNVG3 , ODB , OBP , SPD , 

DIR, LIST 

CALL  INITAIR 
10  CONTINUE 

READ  CONTROL 

IF(AMAXIT.LE.O)  STOP ' END  OF  RUN* 

PRINT  CONTROL 
CALL  AIRMOV 
GO  TO  10 
END 


Subroutine  INITAIR  initializes  most  of  the  variables  for  the  air  flow  simulation. 
In  particular,  it  includes  the  statements  that  read  the  user  data  for  the  zone 
(line  56)  and  opening  (line  74)  descriptions.  Wind  and  stack  pressures  are  set 
to  zero  (lines  91,  92).  Wind  pressure  (PW)  remains  zero  for  all  inter-room 
openings . 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 


C 

C 

C 

c 

c 

c 

c 


SUBROUTINE  INITAIR 

REAL  Z,  T,  F 

INTEGER  MAXAFS 
PARAMETER  (MAXAFS=128) 

INTEGER  MAXZ0N 
PARAMETER  (MAXZON=36) 

COMMON  /AFSL/  FAREA( MAXAFS ) ,FEXP( MAXAFS ) ,ZS( MAXAFS ) ,ZT( MAXAFS ) , 

PW(MAXAFS) , PS (MAXAFS) , FAZM( MAXAFS) , AFSPTR( 2 , MAXAFS) 

INTEGER  AFSPTR 

REAL  FAREA,  FEXP , ZS , PW,  PS,  FAZM,  ZT 

COMMON  /ENVT/  0DB( 1 ) ,OBP( 1 ), SPD( 1 ) ,DIR( 1 ) , AIRDEN, CPAIR, DTR, STDTIM, 
LIST, HCOUNT, ACNVG1 ,ACNVG2 ,ACNVG3 , AMAXIT, FRACT, NAFS , NZON 

INTEGER  STDTIM,  HCOUNT,  AMAXIT,  NAFS,  NZON 
REAL  ODB,  OBP,  SPD,  DIR,  AIRDEN,  CPAIR,  DTR, 
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22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 


ACNVG1 , 


C 


ACNVG2 , 


ACNVG3 , 


FRACT 


COMMON  / ZONL/  MCPM(MAXZON) ,MCPTM(MAXZON) ,TZ(0  jMAXZON) , 

ZZ ( 0 : MAXZON ) , SORTDZ ( 0 : MAXZON ) , PZ(0  : MAXZON) , 
FAHS ( MAXZON ), SUMAF (MAXZON) 


REAL  TZ,  ZZ,  SQRTDZ , PZ,  FAHS,  MCPM,  MCPTM,  SUMAF 
C 

C INITIALIZE  VARIABLES. 

HC0UNT=1 
STDTIM=1 
AMAXIT=20 
FRACT=. 55 
LIST-1 
NAFS=0 
NZON=0 

DTR=0. 0174533 
ACNVG1=. 0001 
ACNVG2=. 000005 
ACNVG3=. 00001 
CPAIR=1004 . 

ODB ( 1 ) =0 . 

0BP( 1 )=100000 . 

SPD( 1 )=5 . 

DIR( 1 )=270. 

SQRT2=SQRT(2.0) 

WRITE(*,101) 

101  F0RMAT( ' 1 ' ,8X, fN ' , 6X, f ZZ ' , 6X, ' TZ ' ,4X, ’FAHS’) 

10  CONTINUE 


C READ  ZONE  DATA: 

C HEIGHT,  TEMP,  SYSTEM  FLOW 


READ  *,  Z , T, F 
IF  (Z.LT.0.0)  GO  TO  20 


NZON-NZON+l 

WRITE(* , 102)  NZON, Z , T, F 

102  FORMATC  ZON: ’ ,I4,3F8.2) 

ZZ(NZ0N)=Z 

TZ (NZ0N)=T 
FAHS ( NZON )=F 
GO  TO  10 
20  CONTINUE 
WRITE(* , 103) 

103  FORMATC 'O' ,8X, *1  N M* ,7X, ' A’ ,7X, ’X' ,7X, 1 C’ ,5X, 'AZM' , 

- 6X, ’ZS’ ,6X, ’ZT’) 

30  CONTINUE 


G 

C 

C 

C 

READ  *,  N , M , A , X , C , AZ 
IF  (N.LE.O)  GO  TO  40 


READ  OPENING  DATA: 

NEAR  SIDE  ZONE,  FAR  SIDE  ZONE, 
AREA,  EXPONENT,  FLOW  COEFFICIENT, 
AZIMUTH,  HEIGHT,  SURFACE  HEIGHT 
Zi  ,Z2 
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73 

NAFS=NAFS+1 

74 

WRITE(*, 104)  NAFS , N , M , A , X , C , AZ , Z 1 ,Z2 

75 

104 

FORMATC  AFS:  ’ , 314 , 3F8 . 4 , 3F8 . 2) 

76 

AFSPTR(1 ,NAFS)=N 

77 

AFSPTR( 2 ,NAFS) =M 

78 

F AREA ( N AFS ) = SQRT2 * C* A 

79 

FEXP(NAFS)=X 

80 

FAZM( NAFS )=AZ 

81 

ZS(NAFS)=Z1 

82 

ZT(NAFS)=Z2 

83 

GO  TO  30 

84 

40 

CONTINUE 

85 

C 

COMPUTE  VARIABLES. 

86 

TZ(0)=ODB( 1 ) 

87 

AIRDEN=0. 0034838*OBP( 1 )/ ( 0DB( 1 )+273 . 15) 

88 

ZZ(0)=0.0 

89 

PZ(0)=0.0 

90 

DO  50  1=1, NAFS 

91 

PS(I)=0.0 

92 

PW(I)=0.0 

93 

50 

CONTINUE 

94 

RETURN 

95 

END 

Subroutine  AIRMOV  computes  the  room  air  densities  (lines  31-33),  the  wind  and 
stack  pressures  for  each  opening  (lines  35-52),  calls  for  solution  of  the  mass 
balance  (line  54),  Computes  the  mass  and  enery  flows  into  each  room  (lines 
56-76),  and  reports  the  room  and  oppening  air  flows  (lines  74-90). 

Descriptions  of  local  variables: 

I index  for  openings. 

N room  number. 

M adjacent  room  number  (0  = ambient). 

V wind  speed  modified  for  height  (m/s). 

W wind  pressure  at  given  height  (Pa). 

Y relative  surface  to  wind  direction  (degrees). 

X wind  pressure  direction  modifier. 

DP  pressure  difference  across  opening  (Pa) . 

MCP  mass  flow  rate  times  specific  heat  of  air  flow  into  room  (W/ C) . 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 


SUBROUTINE  AIRMOV 
REAL  MCP 


INTEGER  MAXAFS 
PARAMETER  (MAXAFS=128) 
C 


c 

c 

c 

c 

c 

c 

c 

c 

10 

c 

c 


INTEGER  MAXZON 
PARAMETER  (MAXZ0N=36) 

COMMON  / AFSL/  FAREA(MAXAFS) ,FEXP (MAXAFS ) ,ZS (MAXAFS ) ,ZT(MAXAFS) , 

PW( MAXAFS) ,PS(MAXAFS) , FAZM( MAXAFS) ,AFSPTR( 2 , MAXAFS) 

INTEGER  AFSPTR 

REAL  FAREA,  FEXP , ZS,  PW,  PS,  FAZM,  ZT 

COMMON  /ENVT/  ODB( 1 ) ,OBP( 1 ) , SPD( 1 ) ,DIR( 1 ) ,AIRDEN, CPAIR,DTR, STDTIM, 
LIST, HCOUNT , ACNVG1 ,ACNVG2 ,ACNVG3 ,AMAXIT, FRACT,NAFS ,NZON 

INTEGER  STDTIM,  HCOUNT,  AMAXIT,  NAFS , NZON 
REAL  ODB , OBP , SPD , DIR,  AIRDEN,  CP AIR,  DTR, 

ACNVG1 , ACNVG2 , ACNVG3 , FRACT 

COMMON  / ZONL/  MCPM(MAXZON),MCPTM(MAXZON),TZ(0:MAXZON), 

ZZ(0 : MAXZON) , SQRTDZ(0 :MAXZON) , PZ(0 : MAXZON) , 

FAHS (MAXZON ), SUMAF ( MAXZON ) 

REAL  TZ,  ZZ,  SORTDZ,  PZ,  FAHS,  MCPM,  MCPTM,  SUMAF 

COMPUTE  ZONE  AIR  DENSITIES. 

DO  10  N=0 ,NZON 

SQRTDZ(N)=SQRT(0.0034838*OBP( STDTIM)/ (TZ(N)+273 . 15) ) 

CONTINUE 

COMPUTE  CONSTANT  DELTA-P. 

DO  20  1=1, NAFS 
N=AFSPTR( 1,1) 

M=AFSPTR( 2,1) 

IF(M.E0 . 0 .AND.  HCOUNT. EQ.l)  THEN 

WIND  PRESSURE. 

V=SPD(STDTIM)*(0.1*ZT(I))**0. 143 
W=0. 5*AIRDEN*V*V 

Y=AMAX1 ( FAZM( I ) ,DIR( STDTIM) ) -AMIN 1 ( FAZM( I ) ,DIR( STDTIM) ) 
IF(Y.GT. 180. ) Y=360. -Y 
IF(Y.LE. 90 . ) X=0. 75-1. 05/90. *Y 
IF(Y.GT. 90. ) X=0. 15/90. *Y-0. 45 
PW(I)=X*W 
END  IF 

STACK  PRESSURE. 

PS(I)=9.80*((ZZ(M)-ZS(I))*SQRTDZ(M)*SQRTDZ(M)- 

(ZZ(N)-ZS(I))*SQRTDZ(N)*SQRTDZ(N)) 

IF( LIST. GE. 1 ) PRINT  * , ' CONSTS : ’ , I ,PS( I) , PW( I) 
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52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 


C 


20 


CONTINUE 


COMPUTE  ZONE  PRESSURES  TO  CONVERGENCE. 


CALL  SOLVZP 

C COMPUTE  ZONE  FLOWS. 

DO  30  N=l,NZON 
MCPM(N)=0.0 
MCPTM(N)=0. 0 
30  CONTINUE 
DO  40  1=1 ,NAFS 
N=AFSPTR( 1,1) 

M=AFSPTR( 2,1) 

DP=PS(I)+PW(I)+PZ(M)-PZ(N) 

IF(DP.GT.O.O)  THEN 

MCP=CPAIR*FAREA( I )* SQRTDZ (M)*DP**FEXP ( I ) 

MCPM(N)=MCPM(N)+MCP 

MCPTM(N)=MCPTM(N)+MCP*TZ(M) 

ELSE  IF(M.GT.O)  THEN 

MCP=CPAIR*  FAREA( I )*SQRTDZ ( N ) * ( -DP ) **  FEXP ( I ) 

MCPM(M)=MCPM(M)+MCP 

MCPTM(M)=MCPTM(M)+MCP*TZ(N) 

END  IF 

40  CONTINUE 
DO  50  N=1 ,NZON 

IF(LIST.GE.l)  PRINT*,  'MCPM:  ' ,N,PZ(N) ,MCPM(N) 

50  CONTINUE 

C IF (LI ST. GE. 1 ) PRINT  FLOWS. 

DO  60  1=1 ,NAFS 
N=AFSPTR( 1,1) 

M=AFSPTR( 2,1) 

F=0.0 

DP=PS(I)+PW(I)+PZ(M)-PZ(N) 

IF(DP.LT.O.O)  THEN 

F=-FAREA( I ) * SQRTDZ ( N ) * ( -DP ) **FEXP ( I ) 

IF(LIST.GE.l)  PRINT*,  'FLOW:  ' , I ,M, FAREA( I ) , SQRTDZ (N) , DP , F 
ELSE  IF(DP.GT.O.O)  THEN 

F=FAREA( I ) *SORTDZ (M)*DP**FEXP ( I ) 

IF(LIST.GE.l)  PRINT*,  ’FLOW:  ’ ,I,M,FAREA(I) , SQRTDZ (M) ,DP,F 
END  IF 

60  CONTINUE 
RETURN 
END 
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Subroutine  SOLVZP  computes  the  room  air  pressures.  It  consists  of  an 
initialization  section  (lines  27-54)  and  an  iteration  section  (lines  56-116). 

The  initialization  is  done  by  computing  the  room  air  pressures  using  flow 
exponents  of  1.  This  produces  a set  of  linear  equations  which  are  solved  by 
Gauss  elimination  (line  50).  In  TARP , initialization  is  done  only  on  the 
first  heat  balance  iteration  of  each  hour  (line  32);  otherwise,  it  is  better 
to  use  the  most  recently  computed  room  pressures.  The  components  of  the 
Jacobian  matrix  (AA)  and  the  net  room  air  flows  are  computed  between  lines  68 
and  95.  The  test  at  line  72  prevents  division  by  zero  at  line  79  and  does  not 
otherwise  effect  the  solution.  The  test  at  line  91  accounts  for  rooms  where 
the  net  air  flow  is  too  small  to  effect  the  room  energy  requirement  and  avoids 
a possible  indefinite  result  (0/0)  at  line  92.  The  primary  convergence  test 
at  line  92  compares  relative  flows.  The  Newton's  method  matrix  is  solved  by 
gauss  elimination  at  line  96.  The  TCOUNT  variable  prevents  endless  iteration 
in  case  of  failure  to  converge  (line  99).  The  simple  improvement  for  successive 
over-corrections  is  in  lines  103  through  106. 


Descriptions  of  local  variables: 


I index  for  openings. 

N room  number . 

M adjacent  room  number  (0  = ambient). 

DP  pressure  difference  across  opening  (Pa) . 

F mass  flow  rate  through  the  opening  (kg/s). 

DF  derivative  of  the  flow  rate. 

AA  Jacobian  matrix. 

BB  sum  of  air  flows  into  room  (kg/s). 

CC  Newton  correction  to  the  room  pressures  (Pa) . 
DD  Prior  iteration  corrections  (Pa). 
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SUBROUTINE  SOLVZP 
INTEGER  MAXAFS 
PARAMETER  (MAXAFS=128) 

INTEGER  MAXZON 
PARAMETER  (MAXZ0N=36) 

COMMON  /AFSL/  FAREA(MAXAFS ) ,FEXP(MAXAFS) ,ZS(MAXAFS) ,ZT(MAXAFS) , 

PW( MAXAFS) , PS(MAXAFS) , FAZM( MAXAFS ) , AFSPTR( 2 , MAXAFS) 

INTEGER  AFSPTR 

REAL  FAREA,  FEXP , ZS,  PW,  PS,  FAZM,  ZT 

COMMON  /ENVT/  ODB ( 1 ) ,0BP( 1 ) , SPD( 1 ) ,DIR( 1 ) , AIRDEN , CPAIR, DTR, STDTIM, 
LIST, H COUNT , ACNVG 1 , ACNVG2 , ACNVG3 , AMAXI T , FRACT , N AFS , NZON 

INTEGER  STDTIM,  HCOUNT,  AMAXI T,  NAFS,  NZON 
REAL  ODB,  OBP , SPD , DIR,  AIRDEN,  CPAIR,  DTR, 

ACNVG1 , ACNVG2,  ACNVG3 , FRACT 

COMMON  / ZONL / MCPM ( MAXZON ) , M CP TM ( MAXZON ) , TZ ( 0 : MAXZ ON ) , 
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45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
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52 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 


ZZ(0:MAXZON) , SQRTDZ(0 rMAXZON) ,PZ(0 :MAXZ0N) , 

FAHS ( MAXZON ) , SUMAF ( MAXZON ) 

C 

REAL  TZ,  ZZ,  SQRTDZ,  PZ,  FAHS,  MCPM,  MCPTM,  SUMAF 
C 

INTEGER  T COUNT 
LOGICAL  CNVG 

REAL  AA(MAXZ0N, MAXZON) , BB( MAXZON ) , CC( MAXZON) ,DD(MAXZ0N) 

IF(LIST.GE. 1)  PRINT  *,  ' INITIALIZATION  ' 

TC0UNT=0 

IF( HCOUNT. EQ. 1 ) GO  TO  30 
DO  10  N=1 ,NZ0N 
DD(N)=1 . 0 
BB(N)=FAHS(N) 

DO  10  M=1 ,NZON 
10  AA(N ,M)=0 . 0 
DO  20  1=1 ,NAFS 
N=AFSPTR( 1,1) 

M=AFSPTR( 2,1) 

AA(N,N)=AA(N,N)-FAREA(I) 

bb(n)=bb(n)-farea(i)*(ps(i)+pw(i) ) 

IF(M.LE.O)  GO  TO  20 
AA(M,M)=AA(M,M)-FAREA(I) 

AA(N,M)=AA(N,M)+FAREA( I) 

AA(M,N)=AA(M,N)+FAREA(I) 

bb(m)=bb(mHfarea(i)*(ps(i)+pw(i)) 

20  CONTINUE 

IF(LIST.GE. 2)  CALL  DUMP AB ( AA,BB ,NZ0N , MAXZON) 

CALL  GAUSSY(AA,BB,CC,NZON, MAXZON) 

DO  25  N=1 ,NZ0N 
PZ(N)=CC(N) 

IF(LIST.GE. 1)  PRINT*,  'PZ:  ',N,PZ(N) 

25  CONTINUE 

C NEWTON  ITERATION. 

30  CONTINUE 

TCOUNT=TCOUNT+l 

IF(LIST.GE. 1 ) PRINT  *,  'BEGIN  ITERATION  ',TCOUNT 
CNVG=. FALSE. 

DO  40  N=1 ,NZON 
BB(N)=FAHS(N) 

SUMAF ( N ) =0 . 0 
DO  40  M=1 ,NZON 
AA(M,N)=0 . 0 
40  CONTINUE 

C EVALUATE  FUNCTIONS  AND 

C PARTIAL  DERIVATIVES. 

DO  50  1=1 , NAFS 
N=AFSPTR( 1,1) 

M=AFSPTR(2,I) 

DP=PZ(M)-PZ(N)+PS(I)+PW(I) 

IF(ABS(DP) .LT.ACNVG3* 
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(ABS(PZ(M) )+ABS (PZ(N) )+ABS(PS( I ) )+ABS(PW(  I ) ) ) ) GO  TO  50 
IF(DP.LT.O.O)  THEN 

F=-FAREA(I)*SQRTDZ(N)*(-DP)**FEXP(I) 

ELSE 

F=FAREA( I ) * SQRTDZ ( M) *DP**  FEXP ( I ) 

END  IF 

DF=F*FEXP(I)/DP 

BB(N)=BB(N)+F 

SUMAF ( N ) = SUMAF ( N ) +AB S ( F ) 

AA(N,N)=AA(N,N)-DF 
IF(M.LE.O)  GO  TO  50 
BB(M)=BB(M)-F 
SUMAF(M)=SUMAF(M)+ABS( F) 

AA(M,M)=AA(M,M)-DF 
AA(N,M)=AA(N,M)+DF 
AA(M,N)=AA(M,N)+DF 
50  CONTINUE 
DO  60  N=1 ,NZON 

IF(ABS(BB(N)) .LE.ACNVG2)  GO  TO  60 
IF(ABS(BB(N)/ SUMAF(N) ) .GT. ACNVG1 ) GO  TO  70 
60  CONTINUE 
CNVG=. TRUE. 

70  CONTINUE 

IF(LIST.GE. 2)  CALL  DUMPAB(AA,BB ,NZ0N ,MAXZON) 

C CHECK  CONVERGENCE. 

IF(CNVG)  GO  TO  999 

IF(TCOUNT.GT.AMAXIT)  STOP ’ ITERATIONS ' 

C SOLVE  AA  * CC  = BB. 

CALL  GAUSSY(AA,BB,CC,NZON,MAXZON) 

C IMPROVE  CONVERGENCE. 

DO  80  N=1 ,NZON 

IF(CC(N)/DD(N).LE.-0.5)  CC(N)=CC(N)*FRACT 
IF( ABS( CC(N) ) .GT.ACNVG3)  DD(N)=CC(N) 

80  CONTINUE 

C REVISE  PZ. 

DO  90  N=1 ,NZON 

PZ(N)=PZ(N)-CC(N) 

IF(LIST.GE.l)  PRINT  *,  ’REVIS:  ' ,N,DD(N) , CC(N) ,PZ(N) ,BB(N) 
, SUMAF(N) 

90  CONTINUE 
GO  TO  30 
C 

999  CONTINUE 

PRINT*,  ’ITERATIONS  ’,TCOUNT 

RETURN 

END 
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Subroutine  GAUSSY  is  a Gauss  elimination  routine  for  computing  the  correction 
terms  in  Newton's  method.  The  elimination  is  done  in  lines  7 through  12. 

There  is  no  pivoting  because  the  matrix  is  diagonally  dominant  and  row  exchanges 
would  never  occur.  Lines  14  through  18  allow  for  an  indeterminate  solution 
which  occurs  if,  for  example,  we  are  solving  for  the  air  flow  between  two 
rooms  but  have  no  openings  from  the  rooms  to  ambient.  The  back  substitution 
is  concluded  in  lines  19  through  23. 

1 SUBROUTINE  GAUSSY(A,B,X,N,MAX) 

2 C 

3 REAL  A( MAX, MAX) , B(MAX) , X(MAX) 

4 C 

5 C 

6 C 

7 

8 
9 

10 
11 

12  10 

13  C 

14 

15 

16 

17 

18 

19 

20 
21 

22  20 

23  30 

24  C 

25  RETURN 

26  END 


GAUSS  ELIMINATION. 
NO  PIVOTING. 

DO  10  K=1 ,N-1 
DO  10  I=K+1,N 
D=A(I,K)/A(K,K) 

B(I)=B(I)-B(K)*D 
DO  10  J=K+1,N 

A(I,J)=A(I,J)-A(K,J)*D 

BACK  SUBSTITUTION. 
IF(ABS(A(N,N)).LT.1.E-12)  THEN 
X(N)=0 .0 
ELSE 

X(N)=B(N)/A(N,N) 

END  IF 

DO  30  I=N-1 , 1,-1 
D=0. 0 

DO  20  J=I+1,N 
D=D+A( I , J)*X( J) 

X(I)=(B(I)-D)/A(I,I) 


Subroutine  DUMPAB  will  dump  the  contents  of  the  Jacobian  matrix  (A)  and  the 
net  room  air  flows  (B)  when  detailed  output  has  been  requested  (LIST=2). 


1 

2 

3 

4 

5 

6 

7 

8 
9 


SUBROUTINE  DUMPAB ( A,B ,N,MAX) 

C 

REAL  A( MAX, MAX) , B(MAX) 

C 

DO  10  1=1, N 

10  WRITE (* , 101 ) I, (A(I , J) ,J=1 ,N) ,B(I) 
101  F0RMAT( 14 , i 1E1 1 . 4 ) 

RETURN 

END 
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APPENDIX  B:  SAMPLE  INPUT  FOR  THE  AIR  FLOW  TEST  PROGRAM 


The  following  sample  input  deck  describes  a nine  room,  single  story  building 
as  shown  in  the  plan  below.  An  opening  is  simulated  on  every  wall.  The  first 
nine  lines  give  the  room  height,  temperature,  and  system  air  flow  (read  at 
line  53  of  INITAIR)  for  each  room.  Room  input  is  terminated  at  the  line  with 
a negative  height.  The  next  24  lines  describe  the  openings  in  the  envelope 
and  partition  walls  (read  at  line  71  of  INITAIR).  Each  line  gives  the  room 
number  (room  numbers  must  be  sequential  from  1),  the  number  of  the  room  on  the 
opposite  side  of  the  wall  (0  = ambient),  the  opening  area,  the  flow  exponent, 
the  flow  coefficient,  the  azimuth  angle  of  the  wall,  the  height  of  the  opening, 
and  the  height  of  the  wall.  Openings  in  partitions  are  described  only  once. 

The  openings  information  is  terminated  by  the  negative  room  number.  The  next 
line  gives  the  convergence,  wind,  and  output  control  data  (NAMELIST  input  read 
at  line  15  of  MAIN).  Simulation  is  terminated  when  the  permitted  air  flow 
iterations  (AMAXIT)  are  zero. 
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AMAXIT=20 ,ACNVG1= .01 ,ACNVG2=. 00001 ,ACNVG3=. 00001 ,SPD=5 .000 ,DIR=0 
AMAXIT=0  / 
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